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ABSTRACT 

The availability of a robust and efficient routine for calculating light curves of a finite source 
magnified due to bending its light by the gravitational field of an intervening binary lens is 
essential for determining the characteristics of planets in such microlensing events, as well 
as for modelling stellar lens binaries and resolving the brightness profile of the source star. 
However, the presence of extended caustics and the fact that the images of the source star 
cannot be determined analytically while their number depends on the source position (relative 
to the lens system), makes such a task difficult in general. Combining the advantages of sev- 
eral earlier approaches, an adaptive contouring algorithm is presented, which only relies on a 
small number of simple rules and operations on the adaptive search grid. By us i ng the para- 
metric representation of critical curves and caustics found bv lErdl & Schneider! (|T993), seed 
solutions to the adaptive grid are found, which ensures that no images or holes are missed. 
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1 INTRODUCTION 

In order to determine the characteristics of planets that reveal their 
existence in microlensing events, to model stellar lens binaries, and 
to distinguish between these configurations, efficient algorithms for 
calculating light curves of extended source stars due to microlens- 
ing of a binary lens system are required. For a single point-li ke lens, 
the m agnification of a point source is an analytic function teinsteinl 
1936), the magnification of a uniformly bri ght circular sourc e can 
be expressed by means of elliptic integrals (Witt & Ma cll994b . and 
is therefore a one-dimensional integral over a semi-analytic func- 
tion in general, which can be approximated by the product of the 
point-source magnification with a characteristic function o f a sin- 
gle v ariable, depending on the source brightness profile (Gould 
1 1994b . For a binary lens, however, the complexity increases enor- 
mously. Whereas a point source affected by a single point-mass 
lens always has two images, with the only exception of lens and 
source being perfectly aligned, binary lenses yield either three or 
five images depending on the source position relative to the lens, 
where three-image and five-image regions are separated by ex- 
tended caustics. This implies that one cannot as easily integrate 
over the source position as in the single-lens case. Moreover, the 
point-source magnification cannot be expressed as analytic or semi- 
analytic function, but needs to be determined by m eans of solving a 
fifth-order complex polynomial JWitt & Maoll 1995b . While a two- 
dimensional integration over the source is possible in principle, one 
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carefully needs to keep track of the caustic and needs to deal with 
a diverge nt integrand as it is approached. While ray-s hooting al- 
gorithms jKavser et afll 19861 : [Schneider & Weifi 1987) overcome 
the need for determining the images by creating a magnification 
map based on mapping the image plane to the corresponding true 
source positions, summing the hits within a source still is a dis- 
cretization of a two - dimen sional integration. On the other hand, 
ISchramm & Kavsed dl987h have found that contour plot routines 
offer a simple way for plotting the images of an extended s ource , 
which avoids the inversion of the lens equation. iDominikl d 1995b 
later showed that in addition to providing plots, the same technique 
can be used to efficiently derive numerical values, and in particular 
to determine the magnification of extended sources as well as of the 
incurred astrometric shift in t he centroid of t heir unresolved imges 
by applying Green's theorem (Dominik 1998]). In fact, determining 
the image area from a contour plot is similar to the ray-shooting 
magnification map approach in mapping an image plane grid to the 
respective source position, but rather than the magnification being 
determined over the source area, an integral in the image plane is 
discretized. 



In this paper, an algorithm is presented that determines the 
magnification from the image contour using an adaptive grid rather 
than a static one. By recognizing the positions of the images of 
the source centre, whi ch can be calculat ed by solving a fifth-order 
complex polynomial dWitt & Maoll 1995b . the fact that holes in the 
images must include the positions of the point-like deflectors, and 
by finding all images stretching over critical curves using a para- 
metric representation based on the classification of the topology of 
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the critical curves and caustics bv lErdl & Schneidedd 19931) . it is en- 
sured that all images for a binary lens are appropriately considered. 

While Sect. E] reviews the basic properties of gravitational 
lensing, Sect. [3] presents the adaptive contouring algorithm and 
Sect. [4] shows how to determine the magnification of an extended 
source star from the image contour line. Sect.|5]discusses the criti- 
cal curves and caustics of a binary point-mass lens in order to derive 
an algorithm for finding a point inside an image stretching over a 
critical curve. Example light curves are shown in Sect.|6]before the 
paper concludes with a short summary in Sect. [7] 



2 GRAVITATIONAL LENSING 

Light received from a source object at distance Ds is bent due to 
the gravitational field of a thin sheet of m atter at a distance Dl with 
surface mass density £(£') by the angle dSchneidenl 19851) 
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For (3 denoting the true source position angle and 6 the apparent 
position angle of its observed images, this implies the lens equation 
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with k(x') — (-Dl#o) 2 T,(Dl8qx'). It provides a surjective map- 
ping of the image position x to the source position y, but the lack of 
injectivity means that a source may have more than just one image. 
At critical points x c , defined by 
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two images merge, so that the number of images changes if and 
only if the source passes a caustic point y c — y(x c ). In general, 
critical points and caustic points form closed curves, known as crit- 
ical curves and caustics unless these degenerate into a point in spe- 
cial cases. The conservation of surface brightness by gravitational 
lensing, i.e. I(x) = I[y(x)], implies that the total magnification 
of a point source is given by 
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so that it diverges if the source comes to lie on a caustic. 

For an extended source, as a consequence of Liouville's theo- 
rem, source and lens contours of same brightness correspond. This 
means that if a source contour is described by an implicit func- 
tion F(y; p) — 0, with p being a parameter vector specifying the 
contour, then all its image contours are given by F[y(x);p] — 0. 
Therefore, a contour plot pr ovides image contour lines without in- 
version of the lens equation 1 Sch ramm & Kav ser 1987). 



3 ADAPTIVE CONTOURING 
3.1 Data structure and algorithm 

A search grid in the image plane needs to be large enough to cover 
all images and dense enough, so that no images or holes in the 
images are missed. While fulfilling these two condition, the grid 
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Figure 1. The 16 different elementary squares, defined by whether each of 
the corners is inside or outside the contour line, where 'inside' is marked 
by a filled black circle and 'outside' is marked by an open circle (filled 
white). While in- or out-squares are considered fully inside or outside the 
contour line and therefore of no further interest, the contour line passes 
through squares of type 'corner in/out' or 'half-in/half-out'. All four edges 
of squares with corners alternating between inside and outside along their 
boundary are cut by the contour line, while it is not clear which edges the 
contour line runs between, making such squares indecisive. 



resolution should be kept as low as possible. Obviously, an adaptive 
grid that just provides higher resolution in regions where this is 
required turns out to be superior to a fixed-resolution grid. 

Such an adaptive grid can be built by hierachically nesting 
squares that represent parts of the image plane. In order to approx- 
imate the position of the contour line, a ray is shot to the corre- 
sponding image position given by the defined mapping, Eq. ((2J, 
and it is noted whether the corresponding true source position falls 
inside or outside the source contour. In this res pect, the proposed 
algorithm resembles the ray-shooting approach jKavser et al.l 19861 
ISchneider&Weil3lll987l) . but the indicator of whether the light ray 
hits the source is kept with the image position. In fact, with the 
conservation of surface brightness, as pointed out in Sect. [2] the 
image position has the same relation (inside or outside) with re- 
spect to the image contour as the corresponding source position to 
the source contour. 

A square with the inside/outside indicators for its corners con- 
stitutes the elementary data structure, out of which all relevant in- 
formation about the image plane is constructed. Being character- 
ized by its four indicators, there are 2 4 = 16 types of elementary 
squares, shown in Fig.[TJ 12 of these elementary squares, namely 
those of type 'corner in/out' or 'half-in/half-out', define part of the 
contour line, which crosses two edges of the square whose corners 
have different status. If two opposite corners are found to be inside 
and the other two corners outside, the contour needs to cross all 
four edges, but it is not clear which edges it connects. Therefore, 
such indecisive squares need to be subdivided in order to improve 
the resolution (this is rule 1 below). 

The length of an edge of a square of depth k is 2~ h a, where 
a denotes a unit size. Squares are nested so that each square Sj,„ 
of depth k either contains 4 subsquares of depth k + 1 each sharing 
a different corner with Smnt or Smn does not have any subsquare. 
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Figure 2. The two operations subdivide and extend. The elements 
that are being added are shown with dotted lines or shaded areas. The size 
of the inside/outside indicators reflects the depth of the respective square. 



The minimal and initial data structure consists of a square of depth 
— 1 with 4 subsquares of depth 0. 

The version of the adaptive contouring algorithm presented 
here uses just two operations on the data structure, subdivide 
and extend, illustrated in Fig. [2] While subdivide creates 4 
sub-squares of depth fc+1 to a specified square of depth fc, extend 
enlarges the range covered by the adaptive search grid by creating 
a new top-level grid with the same centre as the previous one but 
twice the length of an edge. The following conditions require in- 
vocation of the operations extend (condition £) or subdivide 
(conditions 5i to 53): 

• Condition £ : A point at the boundary of the outermost square 
is found to be inside 

• Condition Si: Square is indecisive, i.e. it has alternating in- 
side/outside corners along its boundary 

• Condition S2: Adjoining corners of a square have same status, 
but status changes among evaluated points (corners of subsquares) 
along the connecting edge 

• Condition S3: Square of depth smaller than a specified fc m i n 
is of type 'corner in/out' or 'half-in/half-out' 

The operations following conditions £ , Si and 52 eliminate 
all problematic cases and need to be executed for squares of all 
depths before the contour line can be determined, whereas call- 
ing subdivide on all squares fulfilling condition 53 leads to all 
squares that contain the contour line having a depth of at least fc m i n - 
Application of extend can cause condition £ only to arise in the 
new surrounding square. Similarly, the application of subdivide 
can cause condition Si only to arise in subsquares. In contrast, ad- 
joining squares need to be checked for condition 52, while any of 
the conditions may arise in new subsquares created on the appli- 
cation of extend. For this reason, extend or subdivide are 
executed recursively if conditions £ or Si arise, whereas if modifi- 
cations lead to squares fulfilling conditions 52 or 53, these are put 
on stacks specific to the condition and the depth of the respective 
square. Since squares for any depth need to be checked for condi- 
tion 52, a square fulfilling it needs to be placed on the stack even if 
it has alreadly been placed on a stack related to condition 53. This 
can leave squares on stacks for condition 53 that have already been 
subdivided. If such squares are encountered, they must be ignored 
and removed from the stack. 

In order to find the contour line for a depth fc m i n , from the top 
level to the depth fc m i n — 1, squares are taken from the correspond- 
ing stack related to condition 53 and the operation subdivide is 
executed. This might cause conditions £ and Si to arise and cor- 




Figure 3. Representation of contour at a given resolution depth by set of 
nested squares after the rules of the adaptive contouring algorithm have 
been applied. Filled black circles mark corners that are inside the contour, 
while open (filled white) circles mark those outside the contour. Their size 
reflects the depth of the outermost square this corner belongs to. The squares 
of types 'corner in/out' or 'half-in/half-out' (see Fig.[T} contain the contour 
line and are shaded in grey. 



responding extend or subdivide operations being carried out, 
while squares fulfilling conditions 52 are placed on a respective 
stack. Subsequently, all stacks with squares fulfilling 53 are being 
cleared by pulling the respective squares from them and execut- 
ing the subdivide operation. If the application of extend for 
squares fulfilling condition £ or of subdivide for squares ful- 
filling condition 52 or were carried out, this might have created 
squares fulfilling condition 53, so that the full procedure needs to 
be repeated until this is not the case. After all these operations have 
been applied, the stacks with squares fulfilling condition 53 with 
depth k ^ fc m i n hold the squares from which to determine the con- 
tour line. Fig.[3]illustrates the representation of the image plane by 
the nested squares that results from applying the adaptive contour- 
ing algorithm for a given resolution depth fc m i n - 

As pointed out in detail in Sect. [4] The desired magnification 
of the source star can be obtained from an evaluation of approxi- 
mated contour lines for a series of increasing fc m i n until a desired 
accuracy is reached. Rather than subdiving all squares that consti- 
tute the contour line to a certain given depth, it would be advanta- 
geous to pick those squares that are expected to provide the largest 
gain in improving the estimate of the source magnification. How- 
ever, this would require some advanced bookkeeping. 



3.2 Finding all images and holes 

The essential need to find all segments of the contour line can be 
met by inserting a representation of at least one point inside each 
image and each hole. A point that falls exactly onto a corner is 
represented directly in the data structure. Otherwise, one can use 
the fact that with any point that is not exactly on the contour line, 
there is always a region around it that is either inside or outside. 
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A point can then be represented by the surrounding region which 
can take the form of a single in- or out-square if the point is inside 
a square and by two adjoining in- or out-squares if the point falls 
on their common edge. If 1) for each region enclosed by a contour 
line, at least one point is represented by an inside corner or one or 
two in-squares, 2) at least one point inside each hole (i.e. a region 
surrounded by a contour line that is not enclosed) is represented by 
an outside corner or one or two out-squares and 3) the representing 
squares do not overlap, the adaptive contouring algorithm is guar- 
anteed to find all segments of the contour line, so that no enclosed 
region or hole is missed. 

The condition for the representing squares not to overlap can 
be met by determining a maximal size for the square surrounding 
each point to be inserted. A simple choice is to consider touching 
squares of the same size around this point and each of the other 
points, use the minimal size and round it to the next smaller 2~ k ' a. 
If the point to be inserted into the data structure does not fall into 
the top-level square, one needs to call extend until this require- 
ment is met. Subsequently, starting at the top-level square, while the 
depth is smaller than the minimal depth fcj or not all corners have 
the right status (inside for images, outside for holes), one needs to 
branch into the larger-depth subsquare that contains the point. If no 
subsquares exist, these must be created by a subdivide opera- 
tion. If the point to be inserted falls onto an existing corner, one 
can stop, whereas if the point falls onto an edge, the search needs 
to be split into searches inside a square right and left (for a hori- 
zontal edge) or above and below (for a vertical edge). Such a split 
becomes necessary at most once. 



3.3 Image contours for binary point-mass lenses 

Let us consider the special case of a lens of total mass M, com- 
posed of two point lenses with mass fractions mi and m2 = 1— mi 
that are separated by the angle d #e, where 
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With the choice 9o — (9e, the lens equation, Eq. l[2}, then reads 
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where the origin of the coordinate system coincides with the centre 
of mass of the binary-lens system and the line connecting the two 
components is along the j/i-axis. Such a binary point-mass lens is 
not only a fair model for a stellar binary but also for an isolated star 
orbited by planets, where usually the effect of just a single planet 
dominates. 

If an extended source is completely outside the caustics, its 
image contours surround three images, whereas these surround five 
images if the extended source is completely inside. Every point in- 
side the source has a point image inside one of the images of the 
extended source, so that the images of any point inside the source, 
in particular the source centre, provide a seed for finding the image 
contour lines as pointed out in the previous subsection. These can 
be found straig htforwardly by solving a fifth-order complex poly- 
nomial dWitt& Mad 1995b . If part of the source is inside and part of 
it outside, the number of images differs among its enclosed points, 
so that between one and four images of the extended source can 
occur, where some contour lines cross critical curves, so that the 



enclosed images of the source contain more than a single image of 
points within the source that are inside the caustic. Finding a point 
inside the caustic and determining its images would provide all five 
necessary seeds. While this is not trivial, one can take a slightly 
different approach. As long as the source encloses a cusp, the con- 
tour line encloses at least one image of any point inside the caustic, 
so that the at least three images of the source centre are sufficient 
as seeds. The only remaining case is that the source extends over a 
single fold caustic, so that one of its images stretches over a critical 
curve with no position inside mapping to the source centre. A point 
inside such an image can be found as a point on the critical curve 
that maps to a local minimum of the distance between source centre 
and caustic. Finding such a point inside requires some sophisticated 
analysis of the critical curves and caustics of the underlying binary 
lens. This is pointed out in detail in Sect. [5] Seeds for holes are 
found easily for a binary point-mass lens since these need to in- 
clude one of the positions of the two point-mass constituents. With 
the lens equation, Eq. Q, having poles there, it is sensible to ensure 
that these positions are defined to be outside any closed contour. 

Fig. E] shows the data structure after the seed points have 
been inserted for an example configuration, for which mi = 0.3, 
d — 1.2 and the centre of a circular source of angular radius #e, 
where = 0.2, is located at y = (—0.1,0.45) along with the 
image contour lines, the critical curve, and the source contour with 
the caustic. Around the positions of the two point masses, namely 
(—0.84, 0) and (0.36, 0), marked by crosses, two out-squares rep- 
resenting holes have been inserted (these positions falling onto an 
edge), whereas the three images of the source centre have been sur- 
rounded by in-squares. However, there is a fourth image stretching 
over a critical curve, where a seed has been given by a local mini- 
mum of the distance between the source centre and the caustic, as 
described in detail in Sect. 15.41 



4 IMAGE AREA AND MAGNIFICATION 

The conservation of surface brightness by gravitational lensing im- 
plies that for uniformly bright sources, the magnification equals the 
ratio between the total area of the images and that of the unaffected 
source. This means that the problem of calculating the magnifica- 
tion of a uniformly bright source reduces to that of determining the 

area of its images. 

As exploited bv lDominikl ( ll99ill998l) . Green's theorem al- 
lows to transform an integral over a region R into an integral 
over its boundary dR. Namely, for two functions P(x\,X2) and 
Q{xi,X2) that are continous in R and whose partial derivatives 
dP/dx2 and dQ/dx\ are also continuous in R, 



(8) 



/ (£" £) dxidx2 = J Pdxi + Qdx 2 . 

R dR 

In particular, for the area of R, one has dQ/dxi — 8P/dx2 = 
1, which is fulfilled with the choice P(xi,X2) = —\%2 and 
Q(xi, X2) = \x\, so that the magnification of a uniformly bright 
circular source star with angular radius #e reads 
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2npi 
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where dR denotes the union of the contours of all its images. 

An approximation of the contour line can be obtained from an 
inspection of all squares of types 'corner in/out' or 'half-in/half- 
out'. The intermediate value theorem implies that the continuous 
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Figure 4. (top) Image contours for a circular source star resulting from the 
deflection of light by the gravitational field of a binary lens for an exam- 
ple configuration. The coordinates denote angles on the sky in multiples 
of the angular Einstein radius 9^, defined by Eq. (6), where the coordi- 
nate origin coincides with the centre of mass of the binary-lens system. The 
separation of the two point masses with fractional mass mi = 0.3 (left) 
and 1 — mi =0.7 (right) in units of the angular Einstein radius 6*e has 
been chosen as d = 1.2. The thick line shows the critical curve. The inset 
shows the caustic along with the source contour, where the source centre 
is located at y = (—0.1,0.45) and its angular radius in units of 6^ is 
pi, = 0.2. (bottom) Set of nested squares with inside/outside indicators 
constituting the data structure that represents the image plane after seeds 
for the image and the hole positions have been inserted. Image seeds were 
the three images of the source centre and a point on the critical curve that 
maps to a local minimum of the distance between source centre and caustic. 
Hole seeds were the positions of the two point masses, namely (—0.84, 0) 
and (0.36, 0), marked as crosses. In lightface, the image contour lines are 
indicated. The size of the inside/outside indicators (black for inside, white 
for outside) reflect the depth of the respective outermost square. 



function F[y(x),p] restricted to the edge of a square has at least 
one root if the signs at the comers differ, which means that the 
contour line has to cross the edge. For elementary square of types 
'corner in/out' or 'half-in/half-out', see Fig.Q] there are two such 
edges whose crossings with the image contour line define its en- 
try to and exit from the respective square. Usually contour plot al- 
gorithms determine the cut of the contour line with the edge by 
linear interpolation using the values of F[y(x),p] at the corners. 
However, since the images are dominated by convex contours, this 
choice systematically leads to an underestimation of the enclosed 
area. From practice, it has emerged that simply using the midpoint 
between the corners provides an estimate that on average deviates 
less from the true image area for moderate precisions. It might hap- 
pen that the image contour line runs through squares of different 
depth. For adjacent squares of different depth joining at an edge, 
the contour line has to be put through the edge of the square with 
larger depth, i.e. the smaller one. 

Let us consider a circular source with angular radius p* 8e 
whose center is at y^°\ so that p = (y' ', p*) constitutes the re- 
spective parameter vector. For the square denoted by the running 
index i, let (x^ , x% ) and (x^ +1 \ x^ 1 ^) be the adopted intersec- 
tions of the edges with the contour line. A symmetric discretization 
of Eq. $9^ then yields 



A(y m , 



(4 i) +4 i+1) )(4 i+1) -4 i) ) 



4irpi 



n 



(10) 



The absolute error in the area cannot exceed the sum of the 
area of all squares through which the contour line runs, but it can 
be expected to be much smaller than that. While the ratio e between 
these squares and the approximated enclosed area gives some mea- 
sure for the quality of the approximation, it vastly overestimates the 
typical error. However, a suitable value for e can be chosen from 
testing how its variation affects the final result. 

Forageneral surface brightness I[y (x); j/ ' , p*], one finds 
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where 
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and the brightness of the source, centered at y' ^ reads 
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with I being the average surface brightness and £(p) being the char- 
acteristic profile function depending on the fractional radius p. 



However, I[y(x);y 
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jumps at the source boundary 



p* if £(1) 7^ 0, and its derivative at that posi- 



tio n jumps if (d£/ dp)(l) 7^ 0. Similarly to the definition made 
bv lDominikl dl998h for a linear limb-darkening profile, this can be 
avoided for general monotically falling £(p) by defining 



+e(o)-2^(i) 



I[y(x);y W , Pi< ]=Ix{ 



for which I[y(x);y ( -°\ p 

,(0)| 



for \y(x)-yW\ <p* 

£(o) - C r 



,(14) 



Ih(»)-i/( 'I_ 
for y (a;) — ' | > p* 

> 0. With the value of I[y(x)y (0) , 



for |j/(as) — y w | > p* not affecting the integral, the magnification 
of the source can then be determined as 
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5 CRITICAL CURVES OF BINARY LENSES 
5.1 Topology 

For studying the critical curves of a binary lens, it is favourable to 
center the coordinate system on the component with mass fraction 
mi rather than on the centre of mass as was done in Sect, |3.3l while 
keeping the line connecting the two components is along the j/i- 
axis. This means that for the two point lenses with mass fractions 
mi and mi = 1 — mi that are separated by the angle d#E one 
finds the lens equation 



y(x) = x - mi — - (1 - mi): 



;-(rf,Q) 
-(d,o)p 



(16) 



where the coordinate system has been centered on the primary and 
the line connecting the two components is along the yi-axis. 

For every mass fraction mi, the topology of the critical curves 
and the corresponding caustics u ndergo transitions at ch aracteristic 
separations d c and d w , given by ( lErdl & Schneidejl 19931) 



mi(l 



mil 



1 



cL 



,1/3 



+ (1 - mi) 



1/3 



3/2 



(17) 



This results in the three cases of close binaries (d < d c ), interme- 
diate binaries (d c ^ d ^ d w ), and wide binaries (d > d w ). For an 
equal-mass binary in particular, the trans itions between the topolo - 
gies occur at d c = l/y/2 and d w = 2 dSchneider & Weifilll986h . 
The three different topologies and their transitions are illustrated in 
Fig. [5] where the case mi = 0.3 has been used as example. For 
d — * 00, the critical curves tend towards circles with radii Jm\ 
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Figure 5. The three different topologies for the critical curves (left) and 
caustics (right) of binary point-mass lenses, illustrated for mi = 0.3, along 
with the transitions, which occur at separation parameters d c = 0.717 or 
dw = 1.943. The choice of the coordinate system corresponds to Eq. H6\ . 
Small dots indicate the positions of the two lens objects. 



or VT — mT around each of the individual point mass, whereas the 
corresponding caustics are diamond-shaped with four cusps each, 
where two of these are on the yi-axis, onto which the components 
of the lens binary fall, and tend towards points at the position of 
each point mass. As the two point masses approach, the inner cusps 
on the yi-axis touch, as well as the corresponding critical points, so 
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that the critical curve becomes a lying '8'. For intermediate bina- 
ries, i.e. d c < d < d w , there is a single closed critical curve in 
the shape of the bone, while the corresponding caustic contains 6 
cusps, 2 of them being on the binary axis. Towards smaller sep- 
arations, the fold line bends in and touches itself for d — d c , so 
that one outer and two inner closed curves result for a close bi- 
nary, where d < d c . Correspondingly, the two upper or lower fold 
branches next to the on-axis cusps touch and the caustic detaches 
into a diamond-shaped 4-cusp caustic with two on-axis cusps and 
two triangular 3-cusp caustic with all those cusps being off-axis. As 
d — > 0, the outer critical curve tends to a circle of unit radius around 
the centre of mass, the inner critical curves contract towards a point 
and move towards the centre of mass. The outer critical curve corre- 
sponds to the diamond-shaped caustic, which tends towards a point 
at the centre of mass, whereas the triangular-shaped caustic, corre- 
sponding to the inner critical curves, decrease more strongly in size 
and move towards infinity. 



5.2 Parametric representation 

Using polar coordinates x = rfcos ip, sin u>), lErdl & SchneideJ 



2.5 



( 1993) have pointed out that critical points have to fulfill an equa- 
tion of the form 

Gt2(r; d, mi) cos 2 <p + ai(r; d, mi) cos <p + ao(r; d, mi) = , (18) 

which can readily be solved for cos tp for any given r, so that crit- 
ical curves can be parametrized by the radial coordinate r for each 
branch defined by the different values of cos tp. In order to do so, 
however, one needs to find the ranges for r, for which critical points 
exist. The cos (^-dependence gives two values £ (0, n) or 

tp^ G (n,2n), where tp^ + yr ' = 2n, for —1 < coscp < 1, 
which reflects the symmetry around the axis connecting the two 
lens components. 
With 



G 
H 



r 2 +d 2 
2dr, 

1 _ ml 

2mi(l 



- mil 



(19) 



lErdl & Schneldeil dl993h found 



H2(r; d, mi) 
cii (r; d, mi) 
ao(r; d, mi) 



H 2 g-2d 2 h, 
H(h-2Gg). 
G(Gg-h) + 



(1 



mi l 



(20) 



The existence of critical points requires the discriminant of the 
quadratic equation for cosy?, Eq. d 1 8b , to be non-negative, i.e. 
D = a 2 — 4aod2 ^ as well as —1 ^ cosp ^ 1. With 
Eq. d!8t . one finds the extreme values for cos ip to be assumed for 
P- = a,2 — ai + ao — or P+ = a-2 + a± + ao = 0, respectively. 
This means that the number of solutions for cos tp can only change 
for D = 0, P- = 0, or P+ — 0, so that the respective values for 
r at given (d, mi) define minimal or maximal values of branches 
(r, tp) of the critical curve. For 02 = 0, i.e. r — m\ irrespective 
of d, C = cos ip has the unique solution C = — ao/ai (for 02 = 
and d > 0, one always finds a± 7^ 0) with Eq. J18t . while one finds 

C± — — 0.5[ai ± (a( — 4,aoa2) 1 ^ 2 ]/o2 , 



(21) 
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Figure 6. Structure of the branches of the critical curve (r, ip) indicated by 
the allowed ranges [r m ; n , r max ] for r marked by thick lines, which corre- 
spond to D = 0, P— = 0, or = 0. The line type reflects the matching 
value of cos tp, either or C_ , Eq. i2i\ . While dashed lines indicate a 
limit for r with C+ , dotted lines refer to a limit with C_ , and dash-dotted 
lines with both. The solid part of the D = contour refers to neither. 
Bracket intervals for the critical radii, as listed in Table[T] are indicated by 
thin dotted lines. The lighter shaded region marks the range for which a so- 
lution with C+ exists, whereas the darked shades region marks the range 
for which both solutions with either C+ or C— exist. The classification of 
the branches as CI to C4, II and 12, and Wl to W3 for close, intermediate or 
wide binaries, is shown in Table [2] The transitions between the topologies 
are indicated by thin solid lines while the points where the critical curves 
touch are marked as (d c ,r c ) or (c( w ,r w ). The contours for D = and 
P- = touch at (d t ,r t ). 



Fig. H] shows the allowed ranges for r as a function of d using 
mi = 0.3 as example. However, the structure of the solutions does 
not depend on the adopted mass fraction. While for given (d, mi), 
there is at least one and up to three solutions for D = 0, denoted in 



the following by r^' , r] 2 ' , r^y , and P+ = 0, denoted by r 



(1) J2) 



and r 



(3) 

, there is always just a single solution for P_ = 0, denoted 
by rp_ . The two additional solutions for D — are only present 
for close binaries, whereas the two additional solutions for P+ = 
exist for wide binaries. These critical radii can be obtained by 
interval bisection using the initial bracket intervals listed in TableQ] 
which are also indicated by small dotted lines in Fig. [6] For d = 0, 



one finds rp_ 



1 which reflects the outer critical curve 



tendings towards the constant r 

(2) 
r D 

(3) 1/4 r j 

— m/ for d 



1 (Einstein circle), while r$ = 
rjy' = reflects the inner critical curves contracting at the origin. 
Moreover, r 



D 



where C+ and C_ coincide for D 



4aoa2 = 0. 



0, but this is not a valid value, 
because cos <p is out of the allowed range. 

Table [2] lists the branches of the critical curve defined by the 
allowed range [r m i n , r max ] for r and the value for cos tp, either C+ 
or C_ as defined by Eq. <2U . While CI and C2 correspond to the 
small two inner critical curves, C3 and C4 represent the outer criti- 
cal curve. Branch C4 only exists for d t < d < d c , where d t is the 
separation for which the curves D = and P_ = touch. These 
conditions imply 02 = ao, so that ai = 2ao or a\ = 2a2. For 
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Table 1. Bracket intervals for the critical radii defining the branches of the 
critical curves of a binary point-mass lens. 



Critical radius Bracket interval 



Separation 



r P _ 



i I/ 2 n 
(d, d + 1) 

[d — d w + r w , d) 

i 1/2 i 
(m/ ,r w ] 

(0, min{m|' 2 d, m| }) 

(m\/ 2 d, r c \ 
r !/ 4 i 



d > d„ 
d > d„ 

d < d c 
d < d c 



For given mass fraction mi and separation parameter d, critical radii need 
to fulfill P_ = 0, P-f = 0, or D = 0. These roots can be found by interval 

bisection of the respective bracket interval, and rjj are separated by 

l /2 

the asymptotic solution m,' d for d — > 0. 

Table 2. Branches of the critical curves defined by the range of the radial 
coordinate r and the branch of cos ip, defined by Eq. (21) . 

Branch Range of r cos ip Separation 



CI 


[rg J : 




c+ 


d < d c 


C2 


[rg J : 




C- 


d < d c 


C3 


[r P _ 




c+ 


d ^ d t < d c 


C3 




r Wi 


c+ 


d t < d < d c 


C4 




r P _] 


C- 


d t < d < d c 


11 


[rg\ 


r P+ ] 


c+ 


d c < d $C d w 


12 




r P ] 


C- 


dc S~ d ^ d w 


Wl 




r (l)l 


c+ 


d > d w 


W2 




r (3) l 


c+ 


d > d w 


W3 






C- 


d > d w 



Branches C1-C4 refer to the close-binary case, 11/12 to the intermediate 
binary and W1-W3 to wide binaries, where transitions occur at d c and d w , 
which depend on mi . At dt, the curves D = and P_ = touch, so that 
branch C4 only exists for dt < d < d c . 

d — > 0, C3 tends to the Einstein circle, whereas CI and C2 degener- 
ate into points. For d = d c , the inner critical curves touch the outer 
one at the radius r c = (mid c ) 1//3 , where CI and C3 merge into 
II, and C2 and C4 merge into 12. For the intermediate binary, 12 is 
the branch near the object with mass fraction mi at the coordinate 
origin, whereas II is the far branch. For wide binaries, Wl gives 
the critical curve around the far object, while W2/W3 surround the 
near object, which both tend to circles around the respective object. 
At d = d w , the critical curves described by branches Wl and W2 
merge at radius r w = (mirfw) 1 ' 3 into 12, whereas W3 and 12 are 
identical. 



5.3 Cusps and fold lines 

In order for the determinant of the Jacobian (dy/dx) of the lens 
mapping to vanish, at least one of its eigenvalues must be zero. 
In two dimensions, singularities for which there are no non-zero 
eigenvalues cannot be a generic feature. In fact, such singularities 
(known as umbilics) do not exist for any binary point-mass lens 
dSchneider & Weifilll986l : ferdl & Schneiderll 19931) . Therefore, we 



can restrict ourselves to the case of exactly one eigenvalue being 
zero for dicussing the critical curves of binary lenses. 

Since the critical curve is defined as the set of points where the 
Jacobian determinant of the lens mapping vanishes, Eq. 10, the unit 
normal vector is given as n = D/\D\, where D = V det ( ^ ) . A 
corresponding tangent vector then follows as t x = (— rii, fii), and 
t y = (f^-) t x is a tangent to the caustic. Therefore, the caustic is 
a smooth curve with defined tangent direction (fold singularity) as 
long as t x is not an eigenvector eo of the Jacobian to the eigenvalue 
zero (cusp singularity). A cusp is therefore defined by the scalar 
equation n ■ eo — 0. 

In order to identify cusps as roots of n • eo, one needs to find 
eo as a smooth function along the critical curve. However, based 
on just the local Jacobian 



fdy \ = ( yu 2/12 

\dx) \ 2/12 2/22 



(22) 



where yij — dxi/dxj and j/y = yji, there is no such global solu- 
tion. With 2/HJ/22 — 2/12 = 0, one finds that 



yu = V j/22 = 2/12 = , 



(23) 



i.e. the Jacobian is diagonal if and only if one of the diagonal el- 
ements is zero, which reflects the requirement that one eigenvalue 
has to be zero. Moreover yu + 2/22 = 2, i.e. the trace of the Jaco- 
bian is 2 and the non-zero eigenvalue is 2. Since the condition for 
a critical point implies that 2/112/22 ^ 0, one finds —1 ^ 2/12 ^ 1, 
^ yu < 2, and < 2/22 ^ 2. As long as yi 1 / 0, an eigen- 
vector to the eigenvalue zero is given by e^ 1 ' = (— 2/12, 2/n)> but 
it flips over as yu passes through zero, since 2/12 changes sign at 
that instance. Similarly, e 2 ' = (2/22 , —2/12) is an eigenvector to the 
eigenvalue zero for 2/22 7^ 0, which flips over as 2/22 passes through 
zero. This means that the branches of the critical curves or caustics 
need further subdivision. 

However, the relative positions of the cusps, the points where 
the Jacobian is diagonal and the points where the branches of cos ip 
merge, are unique characteristics for each topology, as shown in 
Fig. [7] Without restriction of generality, one can assume mi ^ 
0.5 (otherwise simply mirror the coordinates with respect to the 
vertical axis). From this figure, one can read off how to determine 
the positions of the cusps as well as those points for which yu = 
or j/22 = by sucessive bracketing using an interval bisection 
algorithm for root finding. 

Choosing e 2 ^ as the adaptation of eo for all A where 2/22 7^ 
0, minimized the required number of subbranches, since several 
positions with 2/22 = coincide with the points where the critical 
curve crosses the a;i-axis or the respective caustic crosses the y\- 
axis. 

5.4 Finding caustic-associated images 

If the source includes a cusp, there is no image area that is discon- 
nected from the image position(s) of the source centre. That case 
can only arise if the source contour crosses a contiguous fold line 
between cusps twice. Since the caustic y c (A) is a smooth func- 
tion of the parameter A between cusps, a circular contour centered 
around y^ can only cross it twice if the distance 

d(A;y<°>)= [y c (X)-y W ] 2 (24) 



V 'j 

exhibits a local minimum, so that 

dd(A;y(°') _ dy c (A) 



dA 



dA 



[2/ c W-y (0) ] =0 



(25) 
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Figure 7. Branches of the critical curves and caustics for the three different topologies (wide, intermediate, close). Perpendicular lines separate the different 
ranges for r and cos ip as listed in Table [2] whereas triangles mark the points for which the Jacobian becomes diagonal, i.e. either j/n =0 (pointing 
downward) or 1^22 = (pointing upward), and cusps are indicated by circles. Filled dots mark the position of the two point masses, where the respective areas 
are proportional to the fractional mass. The nomenclature of the branches corresponds to Table|2] For d ^ dt, branch C4 vanishes, while C3 extends between 
both intersections of the critical curve with the xi-axis or both intersections of the caustic with the j/i-axis, respectively. 



constitutes a necessary condition, which means that the tangent to 
the caustic has to be perpendicular to the line connecting it with 
the source centre. Only such points can therefore constitute the re- 
maining image positions that need to be inserted into the adaptive 
grid in order to ensure that all image contours are found. 

If the parameter A is chosen as length of the critical curve, one 

finds 



_d_ 

dA 



in 



8x2 



n.2 



d 
dxi 



(26) 



with n being is the unit normal vector. 

If the curvature K y (A) of the caustic is a monotonic function 
over the considered range [A m in, A max ], the distance d(X; y' ') can 
only have a single extremum. Therefore, for each caustic segment, 
a potential root of dn y /dA is determined, providing the final sub- 
division. With t y (A) denoting the caustic tangent, the curvature is 
given as 



Ky(X) = 



(27) 



where the prime denotes the derivative with respect to A, as defined 
by Eq. l |26| >. By further differentiation with respect to A, dK a /dA 
arises as analytical function of up to the 4th derivatives of y(x) 
with respect to Xi or X2- 



6 APPLICATIONS AND LIMITATIONS 

Fig. [D shows example light curves for circular sources with differ- 
ent radii and different brightness profiles affected by a binary grav- 
itational lens, characterized by the separation parameter d = 1.2 
and the mass fraction mi = 0.3. For the smaller source radii 
(p* = 0.05 and p* = 0.1), the source transits the caustic and 
moves completely inside, so that separated fold-caustic passages 
produce distinct characteris tic peaks that c an be described by a 
generic profile function (e.g. Dominik 2004). In these cases, the in- 
sertion of a seed image mapping to a local minimum of the distance 
of the source centre from the caustic was essential. In contrast, the 
source never moves completely inside the caustic for = 0.2 or 
p* = 0.5. For pi, — 0.2, there are epochs for which the source 
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crosses the same fold line (as for the smaller radii) or different ad- 
joining fold lines, with or without the cusp in between. Moreover, 
for p* = 0.5 there are epochs for which a various number of up to 
three cusps are enclosed by the source. 

For p* = 0.5, a light curve for a brightness profile correspond- 
ing to maximal limb darkening, i.e. 

?(P) = ^V / W i - (28) 

is shown along with that for a uniformly bright source. With a 
smaller fraction of the total brightness in the outer parts for the 
limb-darkened source, the source magnification shows a smaller 
rise as the source enters or exits the caustic but is larger in between. 

The application of Green's theorem to replace the integration 
over the image area by an integration along its boundary is a very 
efficient approach if the images are moderately distorted, so that 
a large area for the given boundary length is enclosed (a circle is 
optimal). There is significantly less gain, however, for extremely 
strong distortions leading to the enclosed area resembling a line. 
This case is indeed approached if a source star gets very closely 
aligned with a lens star that is only associated with much less mas- 
sive companions such as orbiting planets. Resulting in large peak 
magnifications, such configurations are of some specific interest 
due to their planet-de tection potential dGriest & Safizadeh 1 19981 ; 
iRattenburv e t al. 2002). For modelling the event with the largest 
peak magni fication recorded s o far, OGLE 2004-BLG-343 with 
~ 3000, lDong et al.l j2006l) have derived a more efficient vari- 
ant of the ray-shooting technique. In fact, the current version of 
the adaptive contouring algorithm is significantly slowed down for 
very small impact angles between source and planet-surrounded 
lens star, but, for uniformly bright sources, the computation of a 
single magnification with a relative uncertainty of 5 x 10~ 4 can 
still be carried out in ~ 200 ms on a 600 MFlops machine for mag- 
nifications in the range A ~ 100-1000, while a result is obtained 
about 20 times faster (~ 10 ms) for moderate A ~ 3. It turns out 
that the same choice for the accuracy parameter e (Sect.|4j provides 
a more accurate result for larger magnifications. For limb-darkened 
sources, one needs to add another numerical integration, which 
slows down the computation a lot. Determining the magnification 
of limb-darkened sources by integrating along the image boundary 
therefore becomes inefficient for highly-distorted images. In this 
case, alternative techniques should be used. However, the adaptive 
grid approach can be combined with other evaluations of the im- 
age area or magnification. Moreover, limb darkening usually adds a 
small shift to the magnification, so that all other parameters can be 
well approximated by assuming a uniformly bright source during 
an initial search (allowing faster computation), whereas limb dark- 
ening only needs to be taken into account in a final refinement step. 
For source stars idealized to be uniformly bright, the adaptive con- 
touring algorithm appears to be much faster than ray-shooting for 
moderate magnifications, while providing a smaller gain for huge 
magnifications, which however can still be substantial (~ 4) for 
A ~ 100. Since the computation of ray-shooting magnification 
maps for a fixed (d, q, p*) takes at least some minutes, the evalua- 
tion of a several hundred data points by adaptive contouring usually 
appears to be faster. 



7 SUMMARY AND CONCLUSIONS 

The presented adaptive contouring algorithm provides an efficient 
way for calculating microlensing light curves, i.e. the combined 
magnification of the images as a function of time. While the use of 



an adaptive grid prevents wasting time by densely covering regions 
of little interest, it can be ensured that for binary lenses no image is 
missed. This is achieved by finding the images of the source centr e 
by solving a fifth-order complex polynomial dWitt & Mad 1 19951) . 
realizing that holes in the images must include the position of one 
of the lens components, and by identifying a position inside images 
stretchin g over a critical curve m aking use of its parametrization 
found bv lErdl & SchneideJ dl993h . 

The efficiency of the proposed algorithm can be enhanced by 
using a prioritization scheme for the subdivision of squares con- 
taining the image contour line (double-adaptive contouring) rather 
than using a common resolution depth, by using a better estimate 
on the uncertainty of the calculated magnification, or using a bet- 
ter approximation of the contour line and its enclosed area or the 
source magnification. 

While the performance of the current implementation is im- 
pressive for weakly-distorted images, the computation gets signif- 
icantly slowed down if images are strongly distorted, i.e. nearly 
degenerate into lines. For limb-darkened sources, the evaluation of 
the magnification by an integration along a contour line becomes 
inefficient for such images. However, the adaptive grid approach 
could be combined with other evaluation algorit hms. 

As already pointed out bv lDominikl (1998), the presented ap- 
proach can also be used for calculating the astrometric microlens- 
ing signal, i.e. the positional variation of the centroid of light com- 
posed of the images. 

The identification of cusps and the parametrization of the 
branches of the critical curve in between can also serve as basis 
for algorithms searching for matching binary-lens parameters of 
microlensing features where the source passes over a fold- or cusp 
caustic. 

An implementation in C is available 
from the author and can be downloaded from 

|http : //star-www . st-and. ac ■ uk/ ~md3 5 /Software . html| 
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Figure 8. Example light curves for sources of different radii and with different brightness profiles affected by a binary lens. The components of the lens binary 
have mass fractions mi = 0.3 (left) and 1 — mi =0.7 (right) and their separation in units of the angular Einstein radius 6 E is given by d = 1.2. The source 
is assumed to move uniformly with respect to the lens at an angle of 60° with respect to the binary-lens axis and to pass through = (—0.1, 0.45), where 
the coordinates y refer to the angular source position in units of the angular Einstein radius 8 E with respect to the centre of mass of the lens system. The arrow 
shows the direction of motion, where the source moves by the angle 9 E within the time <b . Circles show the snapshot of the source for different radii as it 
passes through = (—0.1, 0.45), while straight lines limit the region of the plane traced by it and a dotted line shows the trajectory of the source centre. 
A thick line shows the caustic, while the positions of the components of the binary lens are indicated by dots whose sizes reflect the mass fractions. The epoch 
to refers to the closest approach of the source to the centre of mass of the lens binary, whereas the time at which the source reaches y(°> = (—0.1, 0.45) is 
indicated by a vertical dashed line in the plots showing the light curves. The magnification is shown as decrease in magnitude Am = 2.5 log A. 
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